MODULE GRIDVAR
!
! Verwendungszweck:
!
!      Bereitstellung von Hilfsvariablen fuer die Gittergleichung
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   GW_MD      : Gittergewicht: Dichte
!   GW_ME      : Gittergewicht: innere Energie
!   GW_OPA     : Gittergewicht: Opazitaet
!   beta       : Gitterparameter: raeumliche Glaettung
!   tau        : Gitterparameter: zeitliche Glaettung [sec]
!   taurel     : definition von tau relativ zur thermischen zeitskala
!   n          : Punktkonzentration (neuer Zeitpunkt)
!   nA         : Punktkonzentration zum alten Zeitpunkt
!   Res        : Soll-Gitteraufloesung
!   n_tmp      : zeitlich gegelaettete Punktkonzentration n
!
!-----------------------------------------------------------------------

      use config, only : rkind, np

      implicit none

      private
      public :: GRID, GW_MD, GW_ME, GW_OPA, n, nA, n_tmp, Res, beta, tau, taurel


      real(kind=rkind) :: GW_MD, GW_ME, GW_OPA
      real(kind=rkind) :: n(np), nA(np), n_tmp(np), Res(np)
      real(kind=rkind) :: beta, tau(np), taurel

CONTAINS





SUBROUTINE GRID
!
! Verwendungszweck:
!
!     Berechnung von Hilfsgroessen fuer die Gittergleichung
!
!     Die Anwendung der Gittergleichung ist durch zwei Argumente limitiert:
!        Erstens kann man wegen der Verwendung eines logarithmischen Gitters nicht
!        bis auf R = 0 gehen. Wenn man bis R = 0 rechnen moechte dann muss man am inneren
!        Rand zumindest einen - besser jedoch mehrere (10-20 ?) - Lagrange Punkte setzen.
!
!        Zweitens gibt es np Gitterpunkte, jedoch nur np-1 skalare Werte als Gittergewichte.
!        Daher muss man fuer einen Punkt eine andere Gittergleichung als den Gleichverteilungsansatz verwenden
!
!        Theoretisch koennte man auch die Werte aus den Pseudozellen als Gittergewichte verwenden. Die dortigen Werte
!        (durch RB gegeben) waeren jedoch unter Umstaenden sehr kontraproduktiv.
!        (z.b. wenn rho am Rand mit Gradienten gegen 0 geht)
!
!     Um beide Bedingungen auf einmal zu erfuellen wird der innerste Gitterpunkt (4) Lagrange'sch (oder sonst wie) berechnet,
!     im Zweifelsfall kann man dann noch leicht zusaetzliche Punkte dazunehmen...
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   n_glatt_i/im1  : raeumlich geglaettete Punktkonzentration
!   nA_glatt_i/im1 : raeumlich geglaettete Punktkonzentration zum alten Zeitpunkt
!                    Da der Zusammenhang zwischen n_tmp und n_glatt so einfach ist (nA_glatt kann man sich fuers Ableiten wegdenken)
!                    habe ich n_glatt nur noch als lokale Variable definiert - analog im mathematica
!
!-----------------------------------------------------------------------

      use primvar, only : X, XA, MR, MD, ME
      use config,  only : np, rkind, BCflag
      use global,  only : tst
      use geomvar, only : S_vol, S_volA
      use physco,  only : z12, z0, z1, z2
      use matvar,  only : OPAros


      implicit none

      integer          :: idx
      real(kind=rkind) :: n_glatt(np), nA_glatt(np)


!-----------------------------------------------------------------------
!    Hilfsvariablen fuer die Gittergleichung
!-----------------------------------------------------------------------


      !$OMP PARALLEL DO PRIVATE (idx)
      do idx = 2, np-2  ! brauche die n's ueber einen um 1 groesseren bereich da ich noch raemlich glaette

         n(idx)  = z12 * ( X(MR,idx+1)  + X(MR,idx)  ) / ( (z1 + BCflag(idx)) * X(MR,idx+1)  - X(MR,idx)  )
         nA(idx) = z12 * ( XA(MR,idx+1) + XA(MR,idx) ) / ( (z1 + BCflag(idx)) * XA(MR,idx+1) - XA(MR,idx) )

      end do
      !$OMP END PARALLEL DO



      !$OMP PARALLEL DO PRIVATE (idx)
      do idx = 3, np-3

! so definiert werden skalare Gittergewichte einer, um eine halbe Schale weiter innen liegenden Stelle, verwendet.
! am innersten Gitterpunkt geht das nicht mehr -> siehe Header...

         Res(idx) = sqrt( z1 + n(idx)**2 *   ( &
                    GW_MD  *(z12 * (z1/X(MD,idx)     + z1/X(MD,idx-1)  ) * (X(MD,idx)     - X(MD,idx-1)  ) )**2 + &
                    GW_ME  *(z12 * (z1/X(ME,idx)     + z1/X(ME,idx-1)  ) * (X(ME,idx)     - X(ME,idx-1)  ) )**2 + &
                    GW_OPA *(z12 * (z1/OPAros(idx)   + z1/OPAros(idx-1)) * (OPAros(idx)   - OPAros(idx-1)) )**2 )  )


         n_glatt(idx)  = ( (z2 - BCflag(idx-1) - BCflag(idx+1)) * beta + z1) *  n(idx) &
                       - beta * ( (z1-BCflag(idx-1)) *  n(idx-1) + (z1-BCflag(idx+1)) *  n(idx+1) )

         nA_glatt(idx) = ( (z2 - BCflag(idx-1) - BCflag(idx+1)) * beta + z1) * nA(idx) &
                       - beta * ( (z1-BCflag(idx-1)) * nA(idx-1) + (z1-BCflag(idx+1)) * nA(idx+1) )


         n_tmp(idx)    = (z1 + tau(idx)/tst) * n_glatt(idx) - tau(idx)/tst * nA_glatt(idx)


      end do
      !$OMP END PARALLEL DO


END SUBROUTINE GRID

END MODULE GRIDVAR
